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ABSTRACT 

The rise of large-scale structure in the universe depends upon the statistical distribu- 
tion of initial density fluctuations generated by inflation. While the simplest models 
of inflation predict an almost perfectly Gaussian distribution, more-general models 
predict primordial deviations from Gaussianity that observations might yet be sen- 
sitive enough to detect. Recent measurements of the Cosmic Microwave Background 
(CMB) temperature anisotropy bispectrum by the Planck collaboration have signifi- 
cantly tightened observational limits on the level of primordial non-Gaussianity (PNG) 
in the Universe, but they arc still far from the level predicted by the simplest models of 
inflation. Probing levels of PNG below CMB sensitivities will require other methods, 
such as searching for the statistical imprint of PNG on the clustering of galactic halos. 
During the cosmic epoch of reionization (EoR), the first stars and galaxies released 
radiation into the intcrgalactic medium (IGM) that created ionized patches whose 
large-scale geometry and evolution reflected the underlying abundance and large-scale 
clustering of the star-forming galaxies. This statistical connection between ionized 
patches in the IGM and galactic halos suggests that observations of reionization may 
provide another means of constraining PNG. We employ the linear perturbation the- 
ory of reionization and semi-analytic models based on the excursion-set formalism to 
model the effects of PNG on the EoR. We quantify the effects of PNG on the large- 
scale structure of reionization by deriving the ionized density bias, i.e. the ratio of 
the ionized atomic to total matter overdcnsitics in Fourier space, at small wavenum- 
ber. Just as previous studies found that PNG creates a scale-dependent signature 
in the halo bias, so, too, we find a scale-dependent signature in the ionized density 
bias. Our results, which differ significantly from previous attempts in the literature to 
characterize this PNG signature, will be applied elsewhere to predict its observable 
consequences, e.g. in the cosmic 21cm background. 

Key words: cosmology: theory, large-scale structure of the Universe, inflation, reion- 
ization, first stars, galaxies: statistics 



1 INTRODUCTION 

Observational probes of the epoch of reionization (EoR) are 
of great interest not only because they promise to yield new 
information on primordial galaxies, and the effects of their 
radiation backgrounds on the inter-galactic medium (IGM); 
these probes may some day deliver a wealth of cosmological 
information as well. The power spectrum of redshifted 21cm 
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brightness temperature fluctuations from the EoR is a no- 
table example with great potential for cosmological applica- 
tion. Theoretical investigations to date have mainly focused 
on exploiting the separation of the matter power spectrum 
from the influence of astrophysical uncertainties in the EoR, 
which according to linear perturbation theory is made pos- 
sible by anis otropy from peculiar ve locity in the neutral hy- 
drogen gas (jBarkana fc LoebluOOSi ). The possibility of this 
clean separation has stimulated a broad literature on the 
topic of constraining cosmology with future 21cm surveys 
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e. g. see Furlanetto et alj (|200y) and references therein. Also 



see IShapiro et alJ ( 20131 ) for a recent re-examination of this 



technique] . 

One topic currently at the forefront of research in cos- 
mology is the possibility that the initial conditions for struc- 
ture formation were not perfectly Gaussian. While the sim- 
plest inflationary models - the canonical single-field slow- 
roll models - predict an almost perfectly Gaussian distribu- 
tion of initial fluctuations, more general inflationary models 
predict significant deviations from Gaussianity that obser- 
vations might yet be sensitive enough to detect. There is 
therefore great interest in developing new ways to measure 
primordial non- Gaussianity (PNG) since its detection (or 
non-detection) would have profound implications for infla- 
tionary theory. 

Perfectly Gaussian initial conditions can be entirely 
characterized by the power spectrum of initial density fluc- 
tuations, while PNG is characterized by higher-order statis- 
tics. The lowest-order statistic that can distinguish Gaussian 
from non-Gaussian primordial fluctuations is the 3-point 
correlation function, or equivalently, its Fourier transform 
- the bispectrum. In general since the bispectrum is zero for 
the Gaussian case, we can write it as 



B^{ki,k2,k'i) = fKhF{ki,k2,ks) 



(1) 



where F{k\,k2,kz) is a function only of the magnitudes 
of the wavevectors of any three perturbation modes, fci,2,3 
(assuming translational and rotational invariance), /nl is 
a dimensionless parameter which indicates the level of 
departure from Gaussianity, called the "non-linearity pa- 
rameter" JGangui et al. 1994 ;_ Wang fc Kamionkowskill200CI : 



iKomatsu fc Spergell I2001I : iBabich et all |2004 ). and $ is 
the Bardeen potential fluctuation in the matter-dominated 
epoch. The primordial spectra generated by inflationary 
models vary considerably from one model to the next and 
can be quite complicated. Observational efforts have there- 
fore focused on phenomenological templates which cap- 
ture the dominate shapes, or functional dependences on 
wavenumber, of the spectra generated by wide classes of 
inflationary models. 

An important example is the so-called "local" template, 
in which $(r) is obtained from a quadratic transformation 
of the local Gaussian fluctuation field, 4>{r), according tqj 



Mr) 



1 /-local 

+ /nl 



where <> refers to 



average 



over all 



(2) 
space 



llHodges et all Il990l: iKofman et all Il99ll: ISalopek fc Bondl 
19901: 



.Gangui et al.l 11994 I Verde et all 120001 : 

Wang fc Kamionkowskil I2OO0I : iKomatsu fc Spergel l200lh . 
To first order in /^l^', equation ^ gives a primordial 



bispectrum of the following forno 



^ There are two conventions for equation ||2j in the literature: 
1) The CMB convention in which $ is evaluated in the matter- 
dominated epoch, and 2) The large-scale structure convention in 
which <& is linearly extrapolated to the present day. These two 
conventions result in difEcrcnt normalizations of the non-linearity 
parameter, /J^°£=''(LSS) » 1.3/^°^=''(CMB). Here and throughout 
the rest of the paper, we use the CMB convention. 
^ The bispectrum in equation (|3} is more general than equation 



Bi"'^ (fci, k2,ki) 



2f^T\P4ki)P4k2) 



+P^{kl)P^{k3) + P^{k2)P4.[kz)], 



(3) 



where P^{k) the power spectrum of the Gaussian po- 
tential, defined by {(j){k)4>{k')} = {2TifP^{k)SD{k + fc'). 
This local template is important in the phenomenology of 
PNG, since all single-field models of inflation with stan- 
dartfl assumptions predict /J?/;^' = (5/12) (1 - n^) ^ 0.016, 
where Us — 0.96 is the spectral index of t he primor- 
dial power spectrum ([Acauaviva et al. 20031: Maldacenal 



20031: ICreminelli fc ZaldarriagalJool : ISeerv fc LidsevlbOOa 
Chen et al.ll2007l : ICheung et allBoOsI ). Being able to detect 



or limit /nl^' at such a small level would therefore be an im- 
portant test of single-field models. On the other hand, other 
more-general models of inflation exist which predict larger 
/nl'^'i so it is also possible to exclude or restrict such mod- 
els by limiting /nl^' (Henceforth, we shall remove the label 
"local" from /nl and mean /nl^' unless otherwise specifled). 

Observational constraints have been placed on the lo- 
cal template by flnding the range of allowed amplitudes, ex- 
pressed in terms of /nl. For example, CMB anisotropy mea- 
surements by the Wilkinson Microwave Anisotropy Probe 
nine-year data analysis (WMAP9) flnd a 9 5% confidence 
limit of -3 < /nl < 77 (JBennett et al.ll2012l '). Newer CMB 
anisotropy results based upon the first 15 months of data 
from the Planck satellite mission are repo rted to yield /nl = 
2.7±5.8 [68% confidence limit error bars. lAde etal] (|2013l )]. 
This result from the highly sensitive all-sky CMB anisotropy 
experiment of Planck may be close to the limit which can 
be achieved by meas urement of the CMB bispectrum alone 
(|Komatsu fc Spergel 2001; Babich & Zaldarriaga 20(3). To 
confirm these limits and probe /nl further, it is necessary 
to consider other methods. 

Correlation statistics of the CMB temperature 
anisotropy directly probe the initial fluctuations while 
they are still in the linear regime. However, PNG also 
has potentially observable effects on the large-scale struc- 
ture that develops as the initial fluctuations grow to 
the highly non-linear po int of forming halos (see e.g. 
iDesiacgues fc Seliakll2010al . and references therein). For ex- 
ample, since non-Gaussianity to first order adds a skewness 
to the distribution function of filtered initial density fluc- 
tuations, it can significantly impact the abundance of the 
rarest and most massive halos (e.g. iMatarrese et al.l I2OOC 



Verde et al . 2001; Grossi ct al. '2 0071: ILo Verde et al. 



Grossi et al. 2009: Valagcas 2010 : iLoVerde fc Smith 



2008 



2011 



Desiacgues fc Seliakll2010b ), or of low-density cos mic voids 
( Kamionkowski et al.l 120091 ; iD'Amico et al.ll201ll ), both of 



((2J, as it can also be generated in a n umber of models that do 
not involve the latter (see footnote 34 in lKomatsu et al.ll201ll . for 
example) . It is nonetheless customary to refer to this form of the 
bispectrum as the local template. 

^ By "standard" assumptions, we mean that the curvature per- 
turbation is initially in the Bunch-Davies state, and that there 
is no super-horizon evolution of t he curvature pert u rbation from 
a non -at tractor solution. S e e e.g. lAg uUo fc Pa rked ll201lh . iGanq 
1I2OIII). iGanc fc Komatsu l ll20ld ). lAgullo fc Shanderal 



De^_^^aban 



Namioo et al 



"Be 



a)12), D cv et al.l ll2012l). ICarnev et al.l 



20131) . and IChen et al.l ll2013h for models which 



relax these assumptions. 
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which originate from the tails of the density distribution. 
PNG can also leave a strong signature in the halo bias, i.e. 
the ratio in Fourier space of the fractional halo number over- 
density to the fractional matter overdensity, by introducing 
a scale-dependence which originates from coupling between 
large and sm all scales modes i n the initial non-Gaussian 
distribution J Dalal ct alj l2008l : iMatarrese fc Verdel l2008l : 
lAfshordi fc ToUcv 20081 ). For the local template in equation 
((S}, this halo bias has been found to depart significantly 
from the Gaussian expectation by a correction term that 
approaches the form Ab{k) oc /nl(&g — 1)/^^ in the small-fc 
limit, where ba is the expected Gaussian bias (fo r recent 



analytical derivat io ns, see e.g. Desiacques et al. _ 2011bl: 



Smith et al.l I2OI2I: JAdshead et all 120121 : iD'Aloisio et all 
2OI2I : lYokovama fc Matsubaral I2OI2I ). Based upon this 
assumed k~^ scale-dependence of bffc). iGiannantonio et al.l 
((2013) recently constrained /nl to be in the range 
—37 < /nl < 25 (95% limit in their most conservative 
analysis) using both the large-scale clustering of massive 
galaxies and the int egrated S a chs-W ol fe effect [ s ee als o 
Slosar et all (|2008| ). IXia et all (|2010|). IXia et all (|2010D . 



Xia et al 



20 111 ) and IRoss et al.l (|2013l ) for previous 
constraints based on the large-scale clustering of galax;ies] . 

Fluctuations in the density of neutral and ionized hy- 
drogen in the IGM during the EoR created by the energy 
released by the first stars and galaxies presents another op- 
portunity to observe the difference between Gaussian and 
non-Gaussian initial conditions. During this cosmological 
phase, which ended at z > 6, before the universe was a bil- 
lion years old, expanding ionized patches of the IGM were 
created wherever galactic halos formed to fuel their growth, 
so the large-scale structure of this patchiness was correlated 
with that of the halos. The connection between patchiness 
and halo clustering introduces a new kind of bias, the ion- 
ized density bias, i.e. ratio of the fractional overdensity of 
intergalactic H II to the fractional total matter overdensity, 
which can also reflect the difference caused by PNG. The 
theory of this ionized density bias, however, is dependent 
not only on the clustering of halos, but on the further com- 
plications of galaxy formation and radiative transfer which 
determine: (i) How much and what kind of ionizing radi- 
ation is released by galactic halos of different masses, (ii) 
How the radiation ionizes the surrounding IGM, and (iii) 
How that process feeds back on the ability of galactic ha- 
los to form stars and release more ionizing radiation. For 
galaxies to form stars, they must collapse and accrete the 
baryonic component of the IGM along with the dark mat- 
ter that dominates the halo mass, and then they must make 
these baryons gravitationally unstable and self-gravitating 
within the halos, by radiatively cooling the gas below the 
halo virial temperature. A general picture of how this pro- 
cess unfolded in the ACDM universe is the following. 

The first stellar sources of reionization likely formed 
through radiative cooling from coUisional excitation of 
rotational-vibrational energy levels of H2 molecules, within 
halos with masses between M ~ 10^ — 10* M0, and virial 
temperatures Tvir < 10* K - the so-called minihalos. How- 
ever, this early period of star formation {z > 15) is thought 
to have been quenched by a corresponding rise in the UV 
background, since the H2 molecules needed for efficient cool- 
ing in minihalos were easily photo-disociated by UV photons 
in the Lyman- Werner bands of H2 between 11.2 and 13.6 



eV (|Haiman et al.l I2OO0I I . Halos in this mass range would 
also have been suppressed as sources of ionizing radiation if 
they formed in places where the IGM was already ionized. 
The gas pressure of the ionized IGM, which is heated to 
~ 10*K by photoionization, opposes baryonic gravitational 
collapse into minihalo s, a phenomenon so metimes referred 
to as "Jeans filtering" (jShapiro et al.lll994l ) , and pre-existing 
minihalos would also lose their baryonic content to photoe- 
vaporation by the io nizing radiation inside IGM H II regions 
(IShapiro et al.ll2004 ). 

Halos with Tvir > 10* K, on the other hand, could 
radiatively cool their gas through coUisional excitation of 
atomic hydrogen. These are the so-called atomic cooling 
halos (ACHs). This minimum Tvir for ACHs corresponds 
to a minimum halo-mass scale of roughly Mmin ~ 10* M©. 
In fact, like minihalos, even ACHs with masses below the 
Jeans-filtering scale (M ~ 10^ Mq) may have been sus- 
ceptible to negative feedback from IGM photo-heating if 
they formed within already ionized regions. These low-mass 
ACHs ( "LMACHS" ) may also not have been massive enough 
to overcome the IGM pressure forces which act to prevent 
the accr etion of inter-gala ctic gas - the fuel of star formation 
[see e.g. Illiev et al.l ((20071) and references therein] . The pre- 
cise boundary between LMACHs and the high-mass ACHs 
("HMACHs") which are massive enough to be unaffected 
by Jeans filtering is still uncertain. Regardless of the astro- 
physical uncertainties in the theory, reionization was likely 
driven to completion by galaxies in halos that were mas- 
sive enough to be rare and highly biased at the relevant 
redshifts {z > 6). It is therefore natural to expect charac- 
teristic differences between EoR models with Gaussian and 
non-Gaussian initial conditions, and for these differences to 
lead to new observational signatures of PNG. 

There has already been some work on the th e ory o f 
reionization in the context of PNG. ICrociani et al.l (120091 1 
calculated the effects of the modified ACH abundance due 
to PNG on the reionization history and the mean inte- 
grated electron scattering optical depth. Tea, of the IGM. 
iTashiro fc Sugivamal (|201J) investigated how the modified 
ACH abundance can impact the number count of ion- 
ized bubbles observed in future maps of the brightness 
temperature of the redshifted 21cm background from the 
EoR, which they suggest as a potential probe of /nl. 
IChongchitnan fc Silkl (|2012h found that the effects of PNG 
on nonlinear biasing of minihalos significantly increases the 
root-mean-square (RMS) of fluctuations in the 21cm bright- 
ness temperature during the EoR. If minihalos are able to 
main tain their reservoirs of neutra l hydrogen throughout the 
EoR, IChongchitnan fc Silkl (120121 ) claim, their contribution 
to the 21cm RMS may allow detection of /nl ~ 0{1) by 
next-generation radio telescopes like the Square Kilometer 
Array (SKA). The possibility to detect PNG by measuring 
the power spectrum of brightness temperature fluctuations 
in the 21cm backgro und from the EoR was considered by 
Ijoudaki et al] (|201lh (henceforth JDFKS). They reported 
a scale-dependent signature in the 21cm power spectrum 
for PNG which is absent in the Gaussian case. This scale- 
dependent signature is related to that known already for 
the underlying galactic halos, as expected since the latter 
are the sources of reionization. 

JDFKS modeled reionizatio n by semi-numeral simula- 
tions using the SimFast21 code (jSantos et ahluOKj ). modi- 
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fied to take account of local PNG. This me thod builds upon 
the an alytical approximation introduced bv lFurlanetto et al.1 
(|2004|) in which reionization is statistically assumed to oc- 
cur in some region when it has collapsed a large enough 
fraction of its mass to ionize all the at oms in that re- 
:ion. T he collapsed fraction in the model of lFurlanetto et al.1 
200J) is computed analytically by an application of the ex- 
tended Press-Schechter approximation, s ometimes referred 
to as the excursio n-set model (ESM) (jBond et al.l Il99ll : 
iLacev fc Colelll993h . We shall refer to this analytical approx- 
imation, which applies the excursion-set formalism to model 
reionization statistically, as the excursion-set model of reion- 
ization (ESMR). The semi- numerical SimFast21 simulations 
go beyond the statistical approach of the ESMR by creating 
a 3D realization of the initial density fluctuations on a cubic 
mesh in a finite comoving volume prior to reionization, and 
extrapolating them forward in time by linear theory (the 
Zel'dovich approximation), to produce an evolving 3D map 
of the ionization field over time. 

JDFKS defined a statistical quantity they called the 
"bias of ionized regions," which could be used to predict the 
21cm power spectrum in the small-fe limit for a given reion- 
ization model in terms of the power spectrum of the under- 
lying matter density fluctuations. In that limit, the matter- 
density fluctuations are linear and the effect of PNG al- 
lowed by existing constraints is negligible. The JDFKS bias 
quantity is hx{k) = {l/xH)y/Pxx{k)/Ps5{k), where Sh is the 
mean neutral fraction, and Pxxik) and Pss{k) are the power 
spectra of the ionized fraction and matter density fluctuation 
flelds respectively, i.e. {x{k)x{k')) = {2'k)'^ Pxx{k)5D{k + k'), 
and {5{k)S{k')) = {2nfPss{k)5D{k + k'), where So is the 
Dirac Delta function. Using their simulated maps of the 
ionization field for cases with and without local PNG, JD- 
FKS found that bx has a strong fc-dependence in the non- 
Gaussian simulations, while bx is fc-independent in the Gaus- 
sian simulations on large enough scales. Much like prior 
studies of the halo bias, JDFKS found they could fit b^ (k) 
from a given non-Gaussian simulation, as a function of bx 
from the corresponding Gaussian simulation, using a sim- 
ple fitting formula in which [b^ (fc) — bx] scales as k~^ in 
the small-fc limit. Since this scale-dependent correction for 
PNG was found to be proportional to /nl as well, a measure- 
ment of the 21cm power spectrum at small k might enable 
a determination of the value of /nl. Indeed, JDFKS used 
their fitting formula to predict that the SKA and Murchi- 
son Widefield Array (M WA) surveys c o uld de tect /nl ~ 50 
and ~ 100 respectively. iTashiro fc Hoi (|2012l ) subsequently 
used the JDFKS fitting formula again to explore the con- 
straining power of the cross-correlation between CMB tem- 
perature anisotropics and 21cm fluctuations on /nl. 

We are interested here in the general problem of pre- 
dicting the signatures of PNG in the observable proper- 
ties of the EoR. For this purpose, a more fundamental 
quantity is the ionized density bias, the ratio of the ion- 
ized atomic to total matter overdensities in Fourier space, 
6pHii(fc) = 5pjj,j(fc)/5(fc), where 5p^^^{k) is the Fourier trans- 
form of the contrast in the ionized hydrogen mass density, 
5pjjjj = /Ohii/phii — 1. Our goal here is to derive this ion- 
ized density bias, fepHiii for reionization models with PNG 
and relate it to the corresponding quantity for the Gaussian 
case, feJ^Hii- Along the way, we will also derive the related 
bias parameter bx deflned by JDFKS for which they have 



reported the fltting formula described above. As we shall 
discuss, we find significant differences between our derived 
expression and their fitting formula. 

The methodology of our work is as follows: (i) We 
extend the anal ytical excursion- set model of reioniza- 
tion (ESMR) of iFurlanetto et al.l (|2004 ) to include non- 
Gaussian initial conditions with general bispectra. This ex- 
tension is made possible by the derivation of the non- 
Ga ussian collapsed fracti on in iD'Aloisio et al.l (|2012h [see 
also lAdshead et al.l (|2012l ) for a similar but independent cal- 
culation] . (ii) We use our extension of the ESMR to derive, 
for the first time, expressions for the scale-independent and 
-dependent contributions to the non-Gaussian ionized den- 
sity bias, which apply to PNG with general bispectra. (iii) 
We test our derived expressions against more fundamental 
numerical calculations of the ionized density bias using the 
lin ear perturbat i on the ory of reionization (LPTR) developed 
by IZhang et al.l (|2007|), in which the linearized ionization 
rate and radiative transfer equations are solved in the long- 
wavelength limit for the perturbations to the ionized density 
field. This approach is ideally suited to the problem at hand 
since it directly employs the relevant physics equations gov- 
erning the ionization state of the IGM and is e xpected to be 
most accurate on the large scales of interest (|Zhang et al.l 
i2007l) . The LPTR provides a powerful and computationally 
cheaper alternative to fully non-linear cosmological radia- 
tive transfer simulations, which at the present day are pro- 
hibitively expensive in this context, since they would have 
to both resolve scales down to the smallest galactic halos re- 
sponsible for reionization, and be large enough, approaching 
~ 1 Gpc'' in volume, to capture the mode-coupling effects 
of PNG responsible for the scale-dependent halo bias, (iv) 
Finally, we derive the JDFKS bias parameter from first prin- 
ciples using our non-Gaussian extension of the ESMR. The 
significant differences we find between our result and the JD- 
FKS fitting formula, which we also confirm using the LPTR, 
will lead to our revision of their forecasted 21cm power spec- 
trum cons traints on /nl - the topic of a follow-up paper to 
this work JMao et al.il2013l ). 

The remainder of this paper is organized as follows. 
In !j2l we present basic definitions, our non-Gaussian ex- 
tension of the analytical ESMR, and our expressions for 
the scale-independent and -dependent non- Gaussian ionized 
denisty bias. In ^ we summarize the LPTR formalism. In 
ij4] we present numerical results from both the ESMR and 
the LPTR as well as comparisons between these two ap- 
proaches. In ^we derive our ESMR prediction for the JD- 
FKS bx parameter, and discuss our differences from their 
fitting formula. Finally, we offer concluding remarks in iJH] 

In the plots presented below, we use a fiducial ACDM 
cosmology with parameters Q,n = 0.28, JIa = 0.72, SI5 = 
0.046, Ho = WOh km s"^ Mpc"^ (with h = 0.7), n, = 
0.96 and as = 0.82, consistent with WMAP7 constraints 
(jKomatsu et al.l I2OIII). We also emp l oy th e linear matter 
power spectrum of lEisenstein fc Hul (| 19991 ) . Although the 
expressions we derive apply to PNG with general bispec- 
tra, we use the above local template to generate all plots in 
the remainder of this paper. All distances are reported in 
comoving units unless otherwise stated. 
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2 ANALYTICAL MODEL OF THE IONIZED 
DENSITY BIAS 

2.1 Basic definitions 

Let pm{r,z) be the total matter density at location r and 
redsliift z. We define the matter density contrast, 



6{r,z) 



Pm{r,z) 



Pm 



(4) 



wliere pm denotes the mean matter density. 

Similarly, let pmir^z) and piiii{r,z) denote the mass 
densities of neutral and ionized hydrogen respectively, and 
let pii{r, z) = pYn{r, z)+piiii{r, z) denote the total hydrogen 
mass density. We will from here on refer to pmiir, z) as the 
ionized density. We define the ionized density contrast as 



5puii{'r,z) 



p-aii{r,z) 

PHIl(z) 



(5) 



where pYai{z) is the mean ionized density at redshift z. The 
central quantity that we will use in this work to quantify the 
impact of PNG on the large-scale structure of reionization 
is the ionized density bias, 



&pHIl(fc, z) = 



6{k,z) ' 



(6) 



where S{k, z) and Sp-^-^-^ {k, z) are the Fourier transforms of 
the matter and ionized density contrasts respectively. We 
now discuss how the ionized density bias is related to other 
quantities sometimes used in the literature. 

Let Xi(r,z) = pniiir , z) / pnir , z) denote the ionized 
fraction, and let S^ir^z) denote the ionized fraction con- 
trast. 



'5x(t-, 



Xi{r,z) 

Xi{z) 



(7) 



where Xi{z) is the spatially averaged ionized fraction at red- 
shift z. On large scales (A; ^ 1), the ionized density bias is 
related to the ionized fraction bias, 



bx{k,z) 



Sx{k,z) 
5{k,z) 



(8) 



through the definition of the ionized fraction 
smoothed on scale R, Xi{R,r,z) = J d^r'W{\r' — 
r\,R) pmi{r,z)/pHir,z), where W{\r' - r\,R) is a 
spherically symmetric filter with characteristic cut-off 
scale R. If the smoothing scale is large, then the above 
equation for Xi{R,r,z) yields the leading-order relation, 
6xik,z) + 5{k,z) — 5pHuik,z), in the limit of fc ^ 1, 
where we have assumed that the total hydrogen density 
fluctuations faithfully trace dark matter fluctuations on 
large scales. Hence, the ionized density bias is related to 
the ionized fraction bias for small k through 



b4k,z)^bpjiii{k,z)-l. (9) 

Equation ((9]) can be used to transform the results for 6pHii 
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presented in this paper, which are always restricted to the 
small-fc regime, to bx as needed. 

Consider a spherical region with radius R. Let 
fcoii{r,Mniin,R,z) bc the fraction of mass within that re- 
gion contained in halos with masses above some threshold 
Mmin, i.e. the collapsed fraction. At times, we will consider 
the contrast in the collapsed fraction. 



Sfcouir,Mniin,R,z) 



fcoll{r,Mmin,R,z) 
/coll (Mmin, 2 ) 



1, 



(10) 



where /coii(Mmin, z) is the mean collapsed fraction at z, and 
a corresponding bias parameter. 



bccouik, Mmin, z) 



5ccoll{k,Mmin,z) 

S{k,z) 



(11) 



where the dependence on R has been dropped under the 
assumption of fc ^ 2tt/R. 

2.2 Generalizing the excursion-set model of 
reionization (ESMR) to include PNG 

In thi s section we extend the ESMR of iFurlanetto et al.l 
(|2004 ) for Gaussian initial conditions to include PNG with 
general bispectra, and use it to derive expressions for the 
Gaussian and non-Gaussian ionized density bias parame- 
ters. These expressions will provide analytical expectations 
against which to compare the numerical LPTR results pre- 
sented in 

2.2.1 General Principles of the ESMR 

The basic postulate of the ESMR is that the local ionized 
fraction within a spherical volume with radius R is propor- 
tional to the number of ionizing photons produced within 
that volume or, equivalently, the local collapsed fraction of 
mass in luminous sources, 



Xi{Mmin,R,z) = CeSMR ./coll (Mmin , -R, z) . 



(12) 



Here, the ("esmr parameter accounts for the efficiency of 
this mass in releasing ionizing photons into the IGM (note 
that the condition of a fully ionized volume implies /coii ^ 
Cesmr)- Rather than try to compute this parameter from 
first principles, we will treat it as a free parameter. When 
numerically evaluating our results in ^ we will fix the value 
of Cesmr by specifying the electron scattering optical depth 
for the model. For simplicity, we assume that ACHs provide 
the only sources of ionizing radiation throughout reioniza- 
tion (see |6] for a discussion of the limitations of this as- 
sumption). The minimum Tyir > 10* K criterion for ACHs 
roughly corresponds to a minimum halo mass of 



Mmin ~ 1.3 X lO'^ Mq 



lO^K 



3/2 



l + Z 

21 



-3/2 



03 



-1/2 



0.7 



/-^mol 

1.22 



-3/2 



(13) 



where A^moi is th e mean molecular weight of the gas (see e.g. 
iHaiman fc Holdeiii2003i '). 
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2.2.2 Gaussian initial conditions 

Consider a spherical region with initial comoving radius R, 
containing mass M ~ pm4.TvR /3. In models with Gaussian 
initial conditions, the excursion set expression for the frac- 
tion of mass within this comoving volume that will end up 
in ACHs at redshift z is|j 



(|2012h applied the non- Markovian extensionj of t he ex- 
cursion set formalism of iMaggiore fc Riottd (|2010al lblla) to 
calculate the collap sed fraction perturbatively [see also 
lAdshead et alj (|2012l 'l for a simila r but independent calc ula- 
tion]. We adopt equation (33) of lD'Aloisio et alJ (|2012l ) for 
the collapsed fraction of ACHs, 



/coii(A'/nim, -R, z) = erfc 



Sc - 5r{z) 



^2[Sn.in{z) - Sr(z)\ 



(14) 



where 5c ~ 1.686 is the critical density in the spherical col- 
lapse model (in an Einstein-de Sitter universe), 

6r(z) = D(z) I dV W {\r - r'\,R) S{r' , z = 0) (15) 

is the density contrast, smoothed on scale R with the 
coordinate-space top-hat filter function, linearly extrapo- 
lateqj to the epoch z, with variance 



Sr{z) = al = (Si) = I .^Ml{k,z)Pi,{k) 



(16) 



and Sminiz) is the variance of density fluctuations smoothed 
on scale -Rmin = (3Mmin/47rpm)^ • Above, D{z) is the linear 
growth factor normalized to unity at z = 0, and we have 
written Sr in terms of the power spectrum of the primordial 
potential fluctuation, P#, and the factor MR{k,z) 



MR{k,z) 



2 k'^T{k)D{z) 



g{Q)W{k,R), 



(17) 



between the potential and the smoothed density contrast 
in the cosmological Poisson equation in the synchronous- 
comoving gauge. Here, g{Q) = (1 -|- Zi)~^D~^{zi), where Zi 
corresponds to the initial epoch, i.e. the limit of large red- 
shift, is the linear growth factor of the potential normalized 
to unity in the matter dominated epoch, evaluated at the 
present day [g{Q) ~ 0.76 in our fiducial cosmology], and 
T{k) is the matter transfer function normalized to unity on 
large scales. 



2.2.3 Extension to non- Gaussian initial conditions 

PNG complicates the analytical calculation of the collapsed 
fraction considerably because the additional correlations 
quantified in the three-point correlation function (and in 
general higher-order correlation functions) make the ex- 
cursion set random walks non-Markovian. iD'Aloisio et al.l 

* Strictly speaking, equation I I14I I applies in the case of Gaussian 
initial conditions and a sharp fc-spacc filtering function. Under 
these conditions, the stochastic evolution of excursion set trajec- 
tories is Markovian, and fully tractable analytically. The use of the 
coordinate space top-hat filter necessitates more sophisticated an- 
alytical techniques resulting in correction terms to equation I I14I I. 
For simplicity, we neglect these terms here. 

^ In contrast to lD'Aloisio et al.l 1 120121) . we do not adopt the usual 
convention in which the density field is linearly extrapolated all 
the way to the present day. The convention we adopt in the cur- 
rent work, where the density field is linearly extrapolated to the 
epoch of interest, makes the redshift dependence in many of the 
summarized expressions more transparent. 



/coll (.M'min, R, z) — /coll + ^fcoW , (18) 

where the non-Gaussian correction i^j 



A/<ron(M.,i„,i?,z) = 
B 



A 



Sr 



1 



"/coll 



+ 



Sr 



3 \ Smin — Sr Sc — 5r J dSmin 

5c — {5c — 5r) coth 



SI^ScSrW a/G,, 



Sr 



Oumin 



^^- si^Sc-Sr)^^^-^""' 



2(5c — 5c5r) 



coth 



5c5r 



Sr 



COmin 



with 



"/coll 



{5c - 5r) 



dSn,in V^{S,-nin - S'fl)^''^ 



exp 



{5c - 5Rf 



2(<S'niin — Sr) 



(19) 



(20) 



The mode coupling effects of PNG between small and large 
scales are encoded in the functions. 



A = (Cm) - (<5|) + 3 (^-Jmin) - 3 (^H^in) (21) 



S = {Sr) + (<5ii<5mm) - 2 {5j,5min) 



C = ((5fl5min) - {5r), 



(22) 



(23) 



(where 5mm = <5ii=B„i„) containing the three-point correla- 
tion functionqj of the smoothed, linearly extrapolated den- 



^ See also IParaniape fc ShethI ||2012|'): IParaniape et"all ll2012h; 
IMusso fc Paraniapel 1 I2OI2I) : IMusso fc ShethI ||2012|) : IMusso et all 
1I2OI3) for other works on non-Markovian extensions of the 
excursion-set formalism. 

^ Here we use the le ading-order non-Gauss ian correction to the 
collapsed fraction of ID'Aloisio et al.l 12013) . They showed that 
this expression is co nsistent with the previo usly obtained non- 
Gaussian halo bias iJDesiacaues et al.ll2011bl ). which is in good 
agreement with results from N-body simu l ations at low redshifts 
l|DesiacQues et al.ll2011a^ . ID'Aloisio et al.l l|2012l ') also obtained a 
next-to-leading order non-Gaussian correction which alters the 
collapsed fraction by a few percent. For simplicity, we neglect 
this correction term here. 

* It is technically the connected two- and three-point correlation 
functions that appear in the above ESMR expressions. However, 
since (i5) = 0, the connected two- and three-point correlation 
functions are equ al to the full two- and three-point correlation 
functions (see e.g. ICoorav fc Shethll2002l ). We therefore drop the 
distinction here. 
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sity field, 



(Sr^Sr^Srs) = / 



d=*ki d^'ka 



Bg,{ki,k2,k3 



(27r)» (27r)3' 
Mr, (fci, z)MR^{k2, z)MR3(bi, z) 



(24) 



where fcs = \Jk\ + fcj + 2fci ■ h-^. Our non-Gaussian exten- 
sion of the ESMR is achieved by plugging equation psp into 
equation ((T5|. 



2.3 Ionized density bias from the ESMR 

2.3.1 Gaussian initial conditions 

We co mpute the ionized dens ity bias in essentially the same 
way as lD'Aloisio et alj (|2012l ') computed the halo biaqj. The 
first step is to write down the ionized fraction contrast, equa- 
tion (O. Assuming Gaussian initial conditions, we find 



Jcoll 
fG 
J coll, 



(25) 



where the subscript, 0, in the denominator denotes the 
global mean collapsed fraction, i.e. evaluation in the limit 
of (5_R — >■ and 5*^ — >■ 0. We then plug in equation (|14[) and 
Taylor expand 5x in Sr about Sr — in the limit of Sr — >■ 0. 
This procedure yields the leading-order relationship between 
the ionized fraction contrast and the initial matter contrast 
in the unevolved density field - the so-called Lagrangian ion- 
ized fraction bias. In order to relate this to the Eulerian bias 
defined in equation ((8|, which relates the ionized fraction 
contrast to the evolved matter de nsity contrast , we a dopt 
the spherical collapse model as in IMo fc White! (|l99y). In 
this model, the initial mass contained within the unevolved 
(i.e. Lagrangian) region is conserved, so the fraction of ion- 
ized mass in the evolved (i.e. Eulerian) region at z is equal 
to the fraction of the initial mass that will end up ionized 
at z. Moreover, in the large-scale limit, the evolved matter 
contrast is to first-order equal to the linearly extrapolated 
initial contrast. The large-scale Lagrangian and Eulerian bx 
are therefore equivalent in this framework. As a final step, 
we convert b^ to &pHii using equation ^ to obtain 



"pHII 



= 1 



2 ain/gu.o 

5c 91nS'niin ' 



(26) 



Hence, the ESMR predicts that the large-scale ionized den- 
sity bias is scale-independent in models with Gaussian initial 
conditions, with an amplitude given by equation 



2.3.2 Scale-dependent ionized density bias from 
non-Gaussian initial conditions 

Here we shall derive the ionized density bias fepHii for non- 
Gaussian initial conditions. For comparison with the Gaus- 
sian bias, we will write fepHii a^s follows: 



'■JpUll 



k(^) 



kW 



f,-Hii + Abyji„-fA&;;H„. 



(27) 



where the labels "i" and "d" refer to scale-independent and 
scale-dependent corrections respectively. Before proceeding 
we make some key simplifications. Since our goal here is to 
compute the large-scale ionized density bias, where we will 
ultimately take the limit of |5i{| <C 1 and Sr ^ 1, we can 
recast equation (|19|l in the much simplified form. 



A /-NG _ ('^min) — •^iSmin^R) 

^/coii — :; 



(28) 



'-'min *JR 
2 



^/coll 
C^'-'niin 



Sr 



Sc — {Sc — Sr) coth 



Sr.Sh 



Sr 



COmin 



Additionally, we will require the mean non-Gaussian col- 
lapsed fraction, for which Sr ^ and Sr — > 0. In this limit, 
equations (|18|) and (|28|l simplify to 



-fNG _ sG I A rNG 

Jcoll.O — /coll,0 + lAJcollfi 



and 



A -fNG _ (<5min) / S, 
^/coll,0 — o 



1 \ "./coll.O 



»^min ^c 1 OOm\T\ 



(29) 



(30) 



respectively. 

As before, we write down the ionized fraction contrast. 



Jcoll 
fG 
J coll. 



f<3 
Jcoll 



A rNG 

^Jcoll.O 



Wcoll,() 



+ 



A rNG 
^/coll 
fG 
Jcoll. 



(31) 



Like the an alogous expressi o n for the halo bias [see equa- 
tion (57) of iD'Aloisio et al.l (|2012l )]. this equation contains 
both scale- independent and -dependent contributions to the 
ionized density bias. We give technical details of our deriva- 
tion in the appendix and summarize the main results here. 
The first term on the right-hand side of equation H3ip gives 
the Gaussian term from the last section, equation (|26p . The 
terms from equation H28p with ((5j^in) contribute a scale- 
independent correction from PNG to the bias, whereas one 
of the terms with (iJ^in^ij) yields a non-zero scale-dependent 
correction. The scale-independent correction is given by 



^ F or Gaussian ini t ial co n ditions, see also McQuinn et al.l 1120051 ) 
and lAlvarez et al.l | |2006|) . iMcQuinn et al.l l l2005r) calculated the 



bias of the overdensity of H II bubble counts with respect to the 
underlying matter d ensity field in the standard ESMR approach. 
lAlvarez et al.l KOOSl ) calculated both the collapsed fraction bias 
and the ionized fraction bias - the latter with their own model of 
reionization. 



Afe:^„= 



(3) 



;{3) 



OpHII 



'pHII 



30cOmiii Oc 



s\ 



■2 
min 



On 



(32) 



(^min) /'S'min dcuotes the skcwucss of dcu- 



where S^ 

sity fluctuations smoothed on the M-^ 

dependent correction is given by 



scale. The scale- 
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Table 1. LPTR notation conventions 



Afe),"H'ii(fc) = 25c (&: 



uW 



G 
pHII 



-(3) 



■^r/n(fc) 

Mmin(A:)' 



where the form factor is 



7r(3) 



(fc) = 



d'^ki 



X MR{ki,z)MR{q,z)B^{k,ki,q), 



(33) 



(34) 



with q = \/fc2 + fc^ + 2fe ■ fei (as before, we use the short- 



-(3) 



-(3) 



hand notation T^^^^ = J^hLr . and TVImin = A1ii=R,-„i„)- 
The scale-dependence in equation (|33|l is encapsulated in 
the factor J-'^;jj(fc)/7V(min(fc), which applies to PNG with 
general bispectra. We note that a useful simplification of 
equation H33p can be made in the fc < 0.1 Mpc~^ regime of 
the local template, where J'^J^ik) is approximateljo con- 
stant and equal to /nl, and the smoothing kernel, W{k, R), 
in Almin(fc) can be set to unity. In this case, the scale- 
dependent ionized density bias takes the simplified forrn^^l. 



Quantity 


Zhane et al. (20071 


This paper 


Coordinate vector 


X 


r 


Conformal time 


T 


V 


In J/ — In uo 


IJ- 


Xv 


Source emissivity 


s 


3 


Secondary ionization 


K 


/3s.i. 


boost factor 






Mean photon density 


A 


fy 


in units of fin 






Mean source emissivity 


fs 


is 


in units of nn 






In a 


(jj 


V 


Source spectrum 


P 


s 


power-law index 







Ab 



(d) 

>HII 



(fc, z) = 3/nl [&^HIl(^) - l] 



dcilmJ^O 



(35) 



g{0)D{z)k2T{k)' 

which makes explicit the 1/fc^ scaling in the limit of fc <C 1 
for the local template. 



3 LINEAR PERTURBATION THEORY OF 
REIONIZATION (LPTR) 

As we have shown in the last section, the ESMR predicts 
that the large-scale ionized density bias is scale- independent 
in Gaussian models, and acquires both scale-independent 
and -dependent corrections in the case with PNG, in a 
manner analogous to the halo bias. In § 31 we will check 
th ese predictions agai nst a more general method, the LPTR 
of IZhang et al.l (|2007h . which starts with the fundamental 
equations governing the ionization state of the IGM - the 
ionization rate and radiati ve transfer equatio ns. Here we 
summarize the formalism of lZhang et al.l (|2007l ) and our ex- 
tension of it to include PNG, and in the next section present 
numerical results from both the ESMR and LPTR. 



3.1 Fundamental equations 

The main physical quantities which appear in the LPTR 
formalism are as follows ^^ I: 



^^ The mass thresholds corresponding to T^i^ ~ lO^K are around 

j^^min = 10* M0. For the corresponding smoothing scales, 

(3) 
J^ ■ (k) deviates from its asymptotic value of /nl by ^ maxi- 
mum of a few percent at fc ~ 0.1 Mpc~^. 

^^ As wc discuss in JS] below, the ionized density bias in equa- 
tions I I27II , I I32II , and I I35II is not equivalent to the "bias of ionized 
regions" defined by JDFKS, and therefore equation I I35II should 
not be compared directly to their equation (4). In fact, we will 
show in il5]that equation II35I I implies a very different expression 
for that quantity. 

^^ Our not a tion is different from the original notation of 
IZhang et al.l 1 120071) . See Table [T] for the conversions. These con- 
versions are intended to bring the notation of the LPTR equations 



• The total comoving number density of hydrogen atoms 
(ionized and neutral) and of ionized hydrogen, riH ~ 
nH{r,Tj) and nHii = nHii(T',7;) respectively. Here, r is the 
comoving spatial coordinate vector and ry is the conformal 
time. 

• n^. = n.y{r,ri,Xu,^), the comoving number density of 
photons per unit solid angle, d^fJ, around the propagation 
direction ft, per unit frequency parameter, Xi, = In;/ — lui^o, 
where v is the photon frequency, and vq = 13.6eV/(27rft) is 
the ionization threshold of hydrogen. 

• The differential emissivity of the ionizing photon 
sources: the number of photons emmited per unit comoving 
volume, per uni t conformal t i me, p er unit a;^, per unit solid 
angle. Following IZhang et al.l (|2007l l . we write the emissivity 
as j/4:TT, where j = j{r, rj, x^, fJ). 

The local ionization state of the IGM is governed by the 
ionization rate equation, a continuity equation which takes 
into account photoionization and local recombination rates: 



dr] 



+ V ■ {nmiu) 



— {nn — riHii 



qhWhii 



) / dx„ / d^n Uj ),/^ l3s.i.ix^,Xi) 
Jo J a-^iv) 



a^r,) ' 



(36) 



where u is the comoving velocity of a volume-element of 
ionized hydrogen, a{x,^) is the photoionization cross sec- 
tion, as = 2.6 X 10~^''cm^/s is the case-B recombination 
coefficient at an IGM temperature T = lO'^K, and a{rj) is 
the cosmological scale factor. In the photoionization term of 
equation p6|) . the factor l3s.i.{xu,Xi) accounts for secondary 
ionizations due to energetic free electrons produced when 
X-ray photons ionize hydrogen. In what follows, we adopt 
a soft source spectrum which is dominated by UV photons 



into closer correspondence with standard notation in the astro- 
physical literature on radiative transfer and photoionized nebulae. 
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[see text surrounding equation (|4H) ]. such that secondary 
ionizations are negligible, and /3s. i. niay be set to unity. 

The radiation field that drives the photoionization term 
in equation (|36} evolves according to the radiative transfer 
equation, 



drj 



+ n ■ Vn^i 



4n 



cr{x„) 



nmi ) n. 



(37) 



where H{ri) is the Hubble parameter. As in IZhang et al.l 
((20071), we write the main physical quantities in equations 
H36p and (|37p in terms of spatial averages and linear pertur- 
bations in the following way: 

nmi = uh [xtir]) + Ami(T',r7)] 



els represent two limiting cases of source lifetimes. In source- 
model A the lifetimes are assumed to be much shorter than 
the duration of reionization, whereas in source-model B they 
are assumed to be much longer. 

Consider the emissivity j smoothed in coordinate space 
over large scales using a smoothing kernel with characteristic 
scale R (the Fourier space linear perturbation equations in 
i]3.4l are rendered independent of R in the large-scale limit 
where k <^ 2n/R). If we assume that ACHs emit on average 
7 {x,^) ionizing photons per unit a;^, per collapsed hydrogen 
atom, then we may write the smoothed emissivity in source- 
model A as 



ji{r, r], X,,) = "/"^(x,,)— [nnir, r), ^)/coii(M,nin, R, Sr, r?)] . 

(39) 
In contrast, the smoothed emissivity in source- model B can 
be written as 



"H = nn [1 + 5{r,ri)] 

% = uh [^-liv,^'^) + ^-i{''',r],x^,fl)] 



(38) 



Before going into more detail about the equations governing 
the spatial averages and linear perturbations, we describe 
our models of the source emissivity. 



3.2 The source emissivity 

To facilitate comparison with the analytical expectations de- 
rived from the ESMR, we will assume in our LPTR calcu- 
lations that ACHs provide the only sources of ionizing ra- 
diation, so only halos above Mmin given by equation lfT3|) 
contribute to the source emissivity. 

In what follows, we consider two source models. In the 
first model, henceforth referred to as "source-model A," we 
assume that the number of photons produced in some time 
interval Arj is proportional to the change in the number of 
hydrogen atoms in collapsed halos in that time interval. In 
this model, photon-production is fueled only by newly col- 
lapsed hydrogen, which may be a reasonable approximation 
if either internal or external feedback mechanisms quickly 
act to limit continuous star formation, and/or if the steep 
rise in abundance of ACHs results in the effective domi- 
nance of reionization by newly collapsed halos. In the sec- 
ond model, "source-model B," we assume that the rate of 
photon production at a time rj is proportional to the num- 
ber of collapsed hydrogen atoms at that time. In this model, 
photon-production is continuously fueled by hydrogen once 
it collapses into halos. 

Our motivation for considering the two scenarios above 
is that source-model A is more similar to the ansatz adopted 
in the ESMR, where the number of ionizing photons pro- 
duced in some region is assumed proportional to the num- 
ber of collapsed baryons in that region. On the other hand, 
radiative transfer simulations often assume that the pho- 
ton production rate is proportional to the collapse baryon 
number, as in source-model B. It is therefore useful to com- 
pare results using both assumptions. Physically, these mod- 



ifl(r,7?,a;i.) = -y'^ {x„)nH{r,Tj,R)fcoii{Mniin,R,SR,Tj), (40) 

where j^ {x„) is the number of photons emitted per unit 
x,^, per unit ri, p er collapsed hydrogen atom. Following 
IZhang et al.l (|2007h . we parameterize the source spectrum 
with a power law in v, 



Y ' {x„) dxi, 



'Clptr(1 + ^)exp[(s -f l)a;^]da::^, (41) 

where Clptr is the total number of photons per collapsed 
hydrogen atom, per unit a;^, and CifpTR is the total number 
of photons, per unit rj, per collapsed hydrogen atom, per unit 
Xu- The normalization of equation (|41|l has been chosen so 
that L°° "J^'^ix^) da;^ — Clptr under the condition s < —1. 
In our numerical results presented in ^ we fix Clptr ^Y ^^~ 
ing the value of the electron-scattering optical depth. The 
power law index, s, can be used to shift the spectrum to- 
wards the soft (UV) or hard (X-Ray) photons. In this work, 
we restrict ourselves to a soft source spectrum, with s — —3 
[see IZhang et al.l (|2007l ) for LPTR calculations with a hard 
spectrum in models with Gaussian initial conditions]. 

Finally, we note that the LPTR can be sourced by any 
model of the collapsed fraction, including those extracted 
from N-body simulations. However, for the sake of com- 
parison, we use excursion-set equations (|14|l and (|18fl in 
the Gaussian and non-Gaussian models respectively, so the 
statistics of the halo sources in the ESMR and LPTR cal- 
culations are exactly the same. However, we emphasize that 
the LPTR is otherwise independent of the ESMR in how 
it models reionization; the former incorporates large-scale 
physics of radiative transfer, photoionization, and recombi- 
nations in the IGM, while the latter follows from the ansatz 
in equation (|12|l . 



3.3 The equations of the spatial averages 

The spatial averages of equation (|36p and (|37p are 



dxi 
drj 



:47r(l- 



)|da;,^{/3s.i.)e^cW (42) 



;; — Xi<^mi 
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and 



and 



^ = il _|. Ha^ - ^^ (1 
9?7 47r dx,^ a? 



Xi) ^T-G jj 



(2) 



(43) 



dA^ 



OA^ 



respectively, where the quantity Chii = ('t.hii)/("'Hii)^ is the 
clumping factor for hydrogen recombination, and 






{nm){nj){PB.i.) 



and 



a 



(2) _ {nmrij) 
^" ~ {nHi){n-y) 



(44) 



(45) 



are the photoionization clumping factors. In the case of a 
soft source spectrum, in which secondary ionizations are 
negligible, C^jj « C'^. From here on we will drop the dis- 
tinction and denote both photoionization factors with C^h- 
These clumping factors cannot be calculated analytically, 
since they are sensitive to non-linear density fluctuations on 
small scales, and detailed feedback effects of the radiation 
background. We therefore employ a few simple models for 
the clumping factors, motivated by numerical results previ- 
ously reported in the literature, which we describe in ii4.1l 
Note that the source emissivity appears in the above 
spatially averaged equations through the function ^s . In ap- 
pendix [B] we show that ^s is given by 



^A A, -,9/coll,0 



drj 



in source-model A, and 



e 



B I 

7 {x. 



)/c< 



(46) 



(47) 



in source-model B. For Gaussian initial conditions, we in- 
sert equation (|14|l with 5r = Q and Sr = into the above 
equations. For non-Gaussian initial conditions, we use (I29p 
along with equation (pOt . 

3.4 The linear perturbation equations 

Following I Zhang et al.l (|2007l ). we substitute the conformal 
time with the variable y = lna(77), and define ct{ii,Xv) = 
a{xv)nB_/{Ha^) and de = aBnii/{Ha^) for notational con- 
venience. Under these substitutions, a represents the proba- 
bility that a photon propagating through a neutral Universe 
with frequency parameter Xi, directly ionizes a single atom 
within a Hubble time. Similarly, &b is the average number of 
times a proton in a fully ionized universe recombines within 
a Hubble time. We take the Fourier transforms of equation 
H36p and H37p and keep only terms that are first-order in the 
Fourier transforms of S, As, Ahii, and A^ (from here on de- 
noted by S, Aa, Ahii, and A-, respectively). Noting that the 
dark matter density contrast is to first order proportional 
to the growth factor, D{ri), and that the peculiar velocity is 
proportional to the gradient of the gravitational potential, 
we obtain 



dAt 



dy 



= F2S - FiAmi + I dx^iPs.i.) I d^nA7(l - x^)a 

(48) 



-F3A7 + 



+ 



^i-y (^ 



dy dxi, '' ' ATvHa 
where we have defined the following auxiliary functions 



(49) 



Fi =2&BXi 



(50) 



-I- 47r / di^o-^^ 
Jo 



(/?s.i.)-(l-xO 



dPs 



<5^ — ^i 



„ dluD 

-F2 = — ; Xi 



dy 
+ 4tt I dx^a^^ 



(51) 



-I 



iPs.i.) - {1 - Xi)x, 



dPs 



dcj) 



F-i = {l- Xi)G 



ik ■ f2 

Ha ' 



(52) 



The source emissivity appears in the above equations 
through the linear source fiuctuation, As{k,ri). In appendix 
[Bl we show that 



K{k,v) = ^^7-^(^.)|^{^('7)/coii,o[l + &tcon(fc,77)]}. 

(53) 



in source-model A, and 



Af(fc,77) = 5(fe,7?)7^(x,)/eon,o[l + 6fcoii(fc,77)] (54) 

in source-model B, where the collapsed fraction bias, bfcoii, 
is obtained as described in the last paragraph of appendix 

m 



The procedure fo r numerically solving equations (|48|l 
and (|49|l is detailed in lZhang et al.l (12007|). First, equations 
(|42p and (|43p are solved for the global reionization history, 
which serves as input for the linear perturbation equations. 
In summar y, the linear coordin ate transformation in equa- 
tion (12) of IZhang et al.l ((20071) is used to change equation 
(|43|l into a first-order ordinary differential equation. The re- 
sulting equation can be solved for ^-y, for some initial Xi. The 
solution for ^^ is then used to find a new solution of equation 
(|42|l for Xi. This procedure is repeated until ^-^ and Xi con- 
verge, which tends to happen after about ten iterations. A 
similar procedure is then applied to solve for the linear per- 
tu rbations. The coordinate transformation in equation (12) 
of IZhang et al.l ((20071) is again applied to change equation 
(|49|l into a first-order ordinary differential equation. The re- 
sulting equation is solved for A-y, which is then integrated 
over i1 to obtain the monopole perturbation of the radiation 
field, under the assumption that the monopole of As consti- 
tutes the main contribution. The monopole perturbation of 
the radiation field is then used as input for equation (|48|) . 
and the solutions are iterated until convergence is achieved. 



4 RESULTS 

4.1 Illustrative models of reionization 

Here we describe the five reionization models which are eval- 
uated in this section. The parameters of these models are 
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Table 2. Illustrative models of reionization 



Model 



Chii C-/H Emissivity 



Efficiency 
Parameter" 



LPTRl Eq. [55] 1 

LPTR2 2 1 

LPTR3 10 1 

LPTR4 Eq. [55] 1 

ESMR N/A N/A 



Eq.[39] 


Clptr ~ 70.3 


Eq.[39] 


Clptr ~ 54.2 


Eq.[39] 


Clptr ~ 87.5 


Eq.\M 


C^PTR = 11800 


N/A 


Cesmr = 50.2 



' Here, Clptr gi^^^ t^^ number of photons per collapsed 
hydrogen atom, per unit z^, while Clptr gi^^^ ^^^ number 
of photons per collapsed hydrogen atom, per H^ , per unit 



suirimarized in Table [2] We consider four models using the 
LPTR, in which we must specify the recombination and pho- 
toionization clumping factors, in addition to the source pre- 
scription, i.e. source-model A [eq. (|39l) ] or source-model B 
[eq. (|40|l ]. In the first model, LPTRl, we use source-model A 
a nd Chii vari e s with redshift according to the fitting formula 
in lllievet al] (|2007l ). 

Chii = 26.2917 exp (-0.1822 z + 0.003505 z'^) , (55) 

obtained from a high-resolution, small-scale N-body simula- 
tion, which resolved the mass scale of minihalos (for which, 
e.g., Chii = 3.8, 6, and 10 at z = 15, 10, and 6 respectively). 
They computed the clumping factor from the density field of 
dark matter outside of halos; a procedure which assumes as 
a first-order approximation that the ionized density fluctu- 
ations faithfully trace dark matter fluctuations in the IGM. 
However, the above flt does not account for feedback effects, 
such as photoionization heating, which act to suppress the 
clumping by increa sing the Jeans scale in side ionized regions. 

Most recently, iFinlator et al.l (|2012l ') conflrmed the sup- 
pression of the recombin ation clumping fac t or by photoion- 
ization heating (also s ee [Pawlik et al.l l2009l : iMcQuinn et al.l 
l201ll : lShuU et al.ll2012l ). They found that Chii rises from ~ 1 
at z = 10 to ~ 3 at 2 = 6. However, they did not resolve 
the mass scale of minihalos which dominates the small-scale 
structure responsible for the recombination clumping factor. 
Nevertheless, as an opposite extreme, we adopt a constant 
Chii = 2 in LPTR2, which may be more representative of 
the photo-heated IGM. On the other hand, LPTR3, with 
Chii = 10, serves as an upper bracket to the range of recom- 
bination rates considered in this work. The LPTR4 model 
is used to illustrate differences between sourcing photons 
with the differential and cumulative collapsed fraction [see 
text surrounding equations H39|) and (|40|l ]. It has the same 
clumping factors as LPTRl, but uses source-model B. 

The photoionization cl umping factor requi r ed by these 
models is more uncertain. lAubert fc Tevssieil (|2010|) find 
from small-scale reionization simulations that C^h = 1 
should be a reasonable approximation to the actual value, 
which they find to be suppressed relative to Chii (also see 
iKohler et al.ll2007l ). For simplicity, we set C^h ~ 1 for all of 
the LPTR models in Table [2] Finally, we also evaluate the 
ESMR. 
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Figure 1. Global reionization histories in models with Gaussian 
initial conditions (sec Table [2] for model parameters). The efE- 
ciency parameters in the LPTR models and the ESMR are fixed 
by setting the electron scattering optical depth to Tes = 0.08. 



For a given reionization history, the electron scattering 
optical depth is 



Tes = O-T 



ne(z)--dz, 
dz 



(56) 



where ctt is the Thompson cross section, ne{z) = Xi{z)nu{z) 
is the cosmic-average free electron density, and Zrec is the 



redshift of the recombination epoch [zr, 



lO"*). For aU of 



the above models, we set the normalizations (i.e. by tun- 
ing efficiency parameters, either Clptr or Cesmr) by fix- 
ing Tes ~ 0.08 in the case with Gaussian initial conditions. 
When we consider cases with non-zero /nl, the normaliza- 
tion values in the last column of Table[2]stay fixed, while the 
Tes values vary fractionally by at most a couple of percent 
for the illustrative /nl = ±50 models adopted below, as we 
show in the next section. 



4.2 Global reionization histories 

Before presenting our main results for the ionized density 
bias, it is instructive to consider the global reionization 
histories in the models of Table [2] A reionization history 
is obtained in the LPTR by solving for the mean ionized 
fraction in equations (|42|) and (|43|) . In Figure [T] we show 
the mean ionized fraction as a function of redshift for the 
Gaussian LPTR models of Table [2] We also plot the mean 
ionized fraction from the ESMR, which is just given by 
Xi — Cesmr/o!)1i,o- We note that reionization ends between 
z ~ 8 — 9 in all of the models show n, which may be at odds 
which recent observations [see e.g. [Robertson et al.l (|2013f ) 
and references therein]. We stress, however, that these mod- 
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Figure 2. The eflfects of local PNG on the global reionization 
history in LPTRl. Note that the source efficiency parameter is 
fixed to Clptr ~ '70-?' across all cases shown. Top panel: mean 
ionized fraction as a function of redshift. Bottom panel: ratio of 
non-Gaussian to Gaussian models. 
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Figure 3. The effects of local PNG on the electron-scattering 
optical depth, res, for a range of /nl values. We show the ratios 
of non-Gaussian to Gaussian models (tes is fixed to 0.08 in the 
Gaussian models). 



els are employed below for illustrative purposes to show the 
effects of PNG on the large-scale structure of reionization. 

Figure [2] illustrates how PNG affects the reionization 
history through its impact on the halo abundance. Although 
we display only the LPTRl case for clarity, the results using 
other models are quantitatively similar. The top panel com- 
pares the mean ionized fraction in the Gaussian and non- 
Gaussian models - the latter with /nl = ±50 (we choose 
these /nl values for illustrative purposes). The normaliza- 
tion of the source emissivity is Clptr ~ 70.3 in all three 
cases. The bottom panel shows the ratio of the non-Gaussian 
to Gaussian Xi. The case with /nl = 50(— 50) has an ion- 
ized fraction that is ~ 16% higher (lower) than the Gaussian 
case at 2 ~ 15, when the universe is only a few percent ion- 
ized. However, the differences drop to ~ 3(1)% when the 
universe is about 50(95)% ionized at z ~ 10(8.1). Despite 
the differences early on, all three cases finish reionization at 
approximately the same redshift, z ~ 8. The largest effects 
of PNG in our models occur at higher redshifts during the 
early stages of reionization when, due to the imposed Mmin 
for ACHs [equation p3|) ]. the collapsed fraction above Mmin 
is dominated by halos that are rare at those redshifts. For ex- 



ample, M„ 



4(5.9)(10)x 10^ Mq at 2 = 20(15)(10), . 



sponding to a peak height of v = 5c/^/ S-mi-n ~ 4.5(3.6)(2.6). 
Since PNG alters the reionization history for a fixed 
source emissivity, it also alters the optical depth according to 
equation (|56p . Figure |3] shows the effects of PNG on t^s for a 
range of /nl values. For comparison, models with /nl ~ 100 
can boost the optical depth by up t o only ^ 3%, which i s 
somewhat larger than the results of lCrociani et al.l ()2009f l. 
They used an independent analytical model to calculate a 
~ 1% enhancement for local /nl ~ 100. 



4.3 Large-scale ionized density bias 

Figure U shows the ionized density bias as a function of 
wavenumber for the reionization models in Table (2] with 
Gaussian (top) and non-Gaussian (/nl ~ 50, bottom) initial 
condition^ij. We show results at a single redshift, z — 10.2, 
corresponding to Xi ~ 50% for all models. Figure [5] shows 
the redshift evolution of the ionized density bias for a fixed 
wavenumber, k = 0.01 Mpc~^. The steep decline of fepHii to 
unity at lower redshifts marks the end of reionization, after 



^■^ The comoving particle horizon is R r- 
which corresponds to a wavenumber of k 



4900 Mpc at 2 ~ 10, 
- 10-3 Mpc-i. While 
we have adopted the Newtonian approximation in this work, 
we note that general relativistic effects on the source statistics 
may come into play on scales approaching k ~ 10—^ Mpc" 
[for detailed stud ies on gal a xy b ia s in the c o ntext of gen- 
eral r ela tivity, see IYoo et al.N2009l) : IYooI 1I2OIOI): iBartolo et al.l 
ll201lh: IChallinor fc Lewij mm: iBonvin fc Durreil ll201l[) : 
iBaldauf et al.lll201ll '): lBruni et al.N2012t) : |jeong et al.N2012l) ]. On 
the other hand, the mean free path of UV photons during reioniza- 
tion is always much smaller than the horizon, so we do not antici- 
pate additional general relativistic corrections to the LPTR itself 
to be significant. Prom a practical viewpoint, such large scales are 
usually disregarded when forecasting constraints on cosmological 
parameters from the EoR 21cm power spectrum, because fore- 
ground contami nation is most problematic on those scales (e.g. 
iMao et al.ll2008l) . So general relativistic corrections are unlikely 
to influence the application of our results to the EoR 21cm power 
spectrum. 
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Figure 4. The scale-dependence of the ionized density bias. Top 
panel: for models with Gaussian initial conditions. Bottom panel: 
the scale-dependent signature of local PNG with /nl = 50- We 
show results at a fixed rcdshift of 2 = 10.2, which corresponds to 
a global ionized fraction of Si 50% in all models shown. 



which fluctuations in the IGM ionized density field faith- 
fully trace matter fluctuations (see below for an explanation 
of the steepness). Figure [6] shows the ionized density bias for 
a range of /nl values in the LPTRl model. 

Figures |4] and [5] show that the ESMR matches the 
LPTR2 results particularly well throughout reionization. 
This is perhaps not surprising, since LPTR2 is the model in 
Table [2] for which recombinations play the least significant 
role (Chii = 2). The ESMR also does well at matching the 
evolution of bpuii in the other LPTR models at early times. 



Figure 5. The redshift evolution of the ionized density bias for 
a fixed scale of A; = 0.01 Mpc~^. 



While there are, as one might expect, significant differences 
in the evolution of 6pHii between the ESMR and LPTR mod- 
els at later times, two basic predictions of the ESMR are 
corroborated by the LPTR: 1) The Gaussian ionized den- 
sity bias is scale-independent on large scales. 2) Local PNG 
introduces through the clustering of sources a strong scale- 
dependent signature in the ionized density bias. 

By construction, the ESMR fepHii is determined by 
the source bias, so the similarities in shape between the 
ESMR and LPTR results in Figure |4] confirm, at least at 
Xi ~ 50%, that there is a strong connection between the 
scale-dependence of the source and ionized density bias pa- 
rameters. However, we may ask whether this connection is 



© 0000 RAS, MNRAS 000, 000-000 



14 D'Aloisio, Zhang, Shapiro, and Mao 



— I — I — I — I — I I 1 1 1 
LPTRl 



"T 1 1 1 — I I I I 



z=10.2 



fNL=50 




— fNL=25 



fNL=10 



J I I I I I I I I I I I I I I I I 



0.001 



0.01 
k [Mpc-i] 



0.1 



Ionized Fraction 

1 0.8 0.5 0.3 0.1 0.01 




k=10-4 Mpc-' 

k=10-3 Mpc-i 

k=10-2 Mpc-i 

k=10-' Mpc-i 

I I I I I I L 



10 



15 



20 



Figure 6. The ionized density bias in LPTRl at a fixed z = 10.2 
for a range of local /nl values. 



preserved throughout reionization. The answer is not clear 
a priori, since the LPTR tracks the propagation of ionizing 
radiation, which can, especially towards the end of reion- 
ization, act over long distances to smooth out the patchi- 
ness in the ionized density field, and diminish its correspon- 
dence with the local clustering of sources. We answer this 
question by plotting the ionized density bias for a range of 
scales (fc = 10"", 10"^, 10"^ and 0.1 Mpc"^) as a function 
of redshift. The top and bottom panels of Figure [7] show the 
Gaussian and non-Gaussian LPTRl results respectively. Al- 
though we focus here on LPTRl for clarity, we have checked 
that the results are similar for the other LPTR models in 
Table [1 

In the top panel of Figure [71 all four curves are so sim- 
ilar that they are indistinguishable in the plot, indicating 
that the Gaussian 6pHii is to a very good approximation 
scale-independent in the LPTR for k < 10~^ Mpc"'^ effec- 
tively throughout all of reionization. In the bottom panel, 
the different amplitudes of the curves illustrate the scale- 
dependence of &pHii due to PNG. The important point is 
that the curves decline at approximately the same rate un- 
til the end of reionization, indicating a preservation of their 
relative amplitudes. This is especially true for the top two 
curves, corresponding to the largest scales, with k = 10~" 
and 10"'^ Mpc~^. The small change with time in the relative 
amplitude between the bottom two curves, with k = 10~^ 
and 0.1 Mpc^^, results from the scale-independent non- 
Gaussian correction to the bias, which is only noticeable at 
larger fc , where the s c ale-de pendent term is less dominant. 

As IZhang et al.l (|2007l ) point out, the preservation of 
the halo clustering imprint on the ionized density bias can 
be understood in terms of a characteristic scale - the mean 
free path of ionizing photons through the IGM - which for 
a given xv is 
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Figure 7. Evolution of the ionized density bias across a range 
of scales. Top panel: in the LPTRl model with Gaussian initial 
conditions. All curves here are degenerate, indicating that the 
Gaussian ionized density bias remains to a good approximation 
scale-independent until the end of reionization. Bottom panel: for 
local PNG with /nl = 50. The line styles here match those in 
the top panel. The scale-dependent imprint of the non-Gaussian 
source bias on the ionized density bias is preserved until the end 
of reionization in the LPTR. 



\{x^) = a^ /[nno{xi,){l - Xi)]. (57) 

The relevant scale is the (/-averaged mean free path, 
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Figure 8. The mean free path of photons through the IGM av- 
eraged over frequency u [see equations l|57|l and JSS}]. The u- 
averaged mean free path is the characteristic scale above which 
fluctuations in the ionized density field correspond to fluctuations 
in the source distribution. Below this scale, fepHll is suppressed 
towards unity, since ionized density fluctuations trace matter fluc- 
tuations inside of fully ionized regions. 



(A).= 



J dxi, £,jX{x,y) 



(58) 



which corresponds to the average distance that a photon 
travels before photoionizing a neutral atom. Ionizing radi- 
ation is unable to smooth out the patchiness of the ion- 
ized density field on scales much larger than (A),^. On these 
scales, the ionized density field traces large-scale fluctuations 
in the sources, so 6pHii is expected to have approximately 
the same scale-dependence as the source bias (though with 
a different amplitude). On the other hand, for scales less 
than {X)v, but larger than the Jeans scale, the ionized den- 
sity field typically traces the matter density fluctuations, 
which suppresses the ionized density bias towards 6pHii ~ 1- 
As reionization proceeds, (A)^ increases, so the scale below 
which 6pHii is suppressed increases until bpHii = 1 for all 
scales at the end of reionization. 

Figure [5] shows (A)^ in our four LPTR models. For 
clarity, we show only the Gaussian models, since the non- 
Gaussian results are essentially identical. It is now clear why 
the halo clustering imprint is so well preserved in the LPTR 
throughout reionization, and why 6pHii drops so rapidly to 
unity at the end of reionization. The value of (A)^ is well 
below the scales of interest (k < 0.1 Mpc~^ corresponding 
to > 60 Mpc) until the very end of reionization, at which 
point it is a very steeply rising^. In fact, (A)^ increases so 



rapidly at the end of reionization that bpHii drops to unity 
at approximately the same redshift for the four wavenum- 
bers considered in Figure [T] In the next section, we explore 
the relationship between Gaussian and non-Gaussian bpHii 
implied by the source-clustering imprint. 



4.4 Analytical mapping between Gaussian and 

non-Gaussian ionized density bias parameters 

In the last section we saw the ESMR prediction of a strong 
scale-dependent imprint in &pHii from the source bias con- 
firmed. On the other hand, the ESMR does not capture the 
details of the imprint, i.e. the bias amplitude and its evo- 
lution, in all of our LPTR models. Indeed, we should not 
expect it to; after all, the simple ESMR developed in !j2] 
is a one parameter model, whereas each LPTR model con- 
tains four parameters (counting the source spectrum slope, 
s, which we have not varied here). One approach we could 
take is to expand upon the ESMR by re-parameterizing it to 
better match the LPTR results. While this may be a fruit- 
ful topic of future investigation, it is perhaps more useful at 
this stage to test whether the ESMR at least gives the cor- 
rect relationship between the Gaussian and non-Gaussian 
fepHii. To this end, for each LPTR model in Tabled we con- 
sider the accuracy of a "hybrid" calculation, which takes the 
ESMR equations (|32|l and (|33|l for the non-Gaussian fopHii, 
and inserts the numerical values of the Gaussian &pHii from 
the particular LPTR model. We then compare the results of 
this calculation to results from the full LPTR calculation. 

Figure [9] shows the ratios of fopHii obtained from the hy- 
brid and full LPTR calculations. The top panel shows the 
ratios as a function of wavenumber for a fixed z = 10.2, while 
the bottom panel shows them as a function of redshift for a 
fixed k — 0.01 Mpc~^. Equations (|32|l and (|33|l best capture 
the relationship between Gaussian and non-Gaussian fcpHii 
in the LPTR2 model. In that case, the hybrid result is within 
5% of the full result until the end of reionization. Perhaps 
more importantly, the flatness of the dashed curve in the top 
panel of Figure |9] shows that the scale-dependence of fepHii 
is well approximated by the ESMR. For the other LPTR 
models, the hybrid calculations work best at earlier times, 
reproducing the LPTR results to within a few percent for 
z > 12 (corresponding to Xi ~ 10 — 20%) after which the dis- 
crepancies tend to rise. Nonetheless, the hybrid calculations 
are accurate to within < 20% across 10~^ ^ k < 0.1 Mpc~^ 
up until the end of reionization for all LPTR models. We 
highlight the mild mismatch in the scale-dependence of fepHii 
between the hybrid and full LPTR calculations. This mis- 
match is most obvious in LPTR3; the model in which re- 
combinations are most significant. In fact, the functional 
dependence of bpHii on k appears to depend mildly on Chii. 
This effect may prove to be minimal, if the lower values of 
Chii (~ 1 — 3) suggested by recent studies are confirmed. 

Given the success of equations p2|) and (|33} in mapping 
to a reasonable accuracy the Gaussian to non-Gaussian bpHii 
from the range of LPTR models considered here, we hope 
that these equations will work equally well when applied to 
more detailed cosmological radiative transfer simulations. 



^■^ We note that the LPTR neglects absorption of photons by 
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Lyinan-limit systems, which would cap the steep rise of the mean 
free path at the end of reionization (see Figure [Sjl . 
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Figure 9. Testing how accurately tlie ESMR predicts tiic rela- 
tionship between Gaussian and non-Gaussian ionized density bias 
parameters in our LPTR models. In the "hybrid" calculation, for 
a given LPTR model, we take the Gaussian bias, fejfjjn, and insert 
it into the ESMR expressions [equations II32I I and II33I I] to obtain 
the non-Gaussian, b^^^^. These results are then compared to the 
full LPTR, calculation, which provides the exact fe^gjj- Above, 
wc show the ratio of hybrid to full LPTR fc^^jj as a function 
of k for fixed z = 10.2 (top), and as a function of z for fixed 
A; = 0.01 Mpc-i (bottom). 



5 COMPARISON WITH PREVIOUS WORK 

In this section we compare our results to the results of JD- 
FKS. We begin by summarizing important det ails of the JD- 
FKS SimFast21 reionization simulations [see ISantos et al.l 
(J201CI ) for details on the SimFast21 code methodology]. 
Their initial density fields, both Gaussian and non-Gaussian, 
were generated on a regular cubic mesh with A'^ — 2048'^ 
cells in a (3 Gpc)'^ volume, corresponding to a halo mass- 
resolution of M « 5 X 10^^ M0. In order to account for 
mass in smaller halos, down to their ACH mass-threshold 
of Mniin ~ IO^Mq, JDFKS employed an analytical sub-grid 
model in which equation (|14p - the excursion-set expression 
for the Gaussian collapsed fraction - is use d to fill in the 
unres olved collapsed mass in each mesh cell (jJoudaki et al.l 
[201l|). At z = 7.5, the redshift at which their analysis was 
based, the IO^Mq threshold corresponds to a peak height 
of !/ = 5c/y/S « 2, while the 5 x 10^^ Mq mass-resolution 
limit corresponds to i^ ~ 4.7. Their reionization simulations 
were therefore sourced mainly by unresolved sub-grid halos, 
since these were by far more abundant than halos resolved 
directly by their mesh. 

As we point out in § [ij JDFKS did not directly cal- 
culate the ionized density bias. Rather, they calculated a 
quantity they called the "bias of ionized regions," whose 
definition has caused some confusion in the literature. They 
defined this bias parameter to be bx = \f^ZIT^7s (Note 
that we use the "hat" notation to distinguish this quan- 
tity from the distinct but related ionized fraction bias, h^. 
See below for the relationship between 6i and hx). Here, 
3^f,f, = f^x^Pss and ^^x = f^P^x, where ft is related 
to the spatially averaged 21cm brightness temperature (Tt) 
by Tb — Tb/xu under the usual assumption that the spin 
temperature is much greater than the CMB temperature, as 
expected during reionization. We emphasize the important 
distinction between bx, used in JDFKS, and the ionized frac- 
tion bias, bx, defined in equation (|8]). These quantities are 
related as follows: bx = \/3^xx/^s5 = ii-/xii)\/Pxx/Pss = 



{x^/xu)^/PsJ7/Psc 



y^3^xx/^SS 
which yields 



■^bx 
xh 



(59) 



Using equation (|9]), we also find a simple relationship be- 
tween bx and the ionized density bias, &pHii, which applies 
on large scales: 



bx{k,z) 



Xi{z) 

xii{z) 



[bpmi{k,z) 



(60) 



JDFKS report that the bx computed from their non- 
Gaussian simulations are to a good approximation related 
to the bx from corresponding Gaussian simulations by the 
following proposed fitting functioio: 



Equations (|32|l and H33p should be similarly tested against 
full simulation results once large enough boxes become avail- 
able. 



^^ Here we restore a factor of g{0) in the denominator of the 
scale-dependent term. JDFKS did not normalize D{0) = 1, so 
g(0) does not appear in equation (4) of their paper; it is absorbed 
into their definition of D{z) ijjoudaki 2013). In contrast, we choose 
to normalize D{0) = 1, so g{0) must appear in the denominator 
of our equation II61I I . 
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Figure 10. Testing the accuracy of our ESMR prediction [eq. 
I|63| l] versus the JDFKS fitting formula [eq. l|6H l] in predict- 
ing the relationship between Gaussian and non-Gaussian bx = 
y/^xx/^SS (The JDFKS bias parameter) in our LPTR calcu- 
lations. To obtain the curve labeled "hybrid", we take fe^ from 
LPTRl and insert it into equation 11631 . The curve labeled "JD- 
FKS" is obtained in a similar way, but using l|61|l . For clarity, we 
only show results from the LPTRl model in Table [2] 



bx (k, z) = bx {z) -I- 3/nl 



b'^{z) 



Ot;! Lfi 



g{Q)D{z)k'^T{k)' 



(61) 

The quantity 5x is analogous to the critical overdensity for 
collapse which appears in the well-known non- Gaussian halo 
bias formula (Matarr ese fc Verdel l2008l : iDalal et all l2008l : 
lAfshordi fc ToUcv 20081 '). The condition fcoii{Mmin, R, z) > 



Cesmr ^'^^ ^ fully ionized region in the ESMR along with 
equation (|14p implies (for Gaussian initial conditions) a min- 
imum linearly extrapolated overdensity of 



5x ^ 5c — V2iC'((^ESMR)\/5'min — Sr 



(62) 



for a region with Lagrangian radius R to self-ionize, where 
A'(Cesmr) = erf-i(l - Cesmr)- JDFKS defined Sx to be 
the "average" Sx- For illustrative purposes, they calculated 
5x = 1.1 by averaging equation (|62|l over the fraction of 
volume filled by H II bubbles. However, they set Sx = 1 for 
convenience when they used equation (|61|l to forecast 21cm 
power spectrum constraints on /nl- 

We note that equation (|61|l was not derived from an- 
alytical theory. Rather, it was postulated by JDFKS based 
on their intuition from previous literature on the halo bias. 
Here, we derive a formula for b^ from first principles in the 
ESMR, and compare it directly to equation (|6ip . We then 
compare the accuracies of our result and equation (I61II in 



capturing the relationship between b^ and bx when both 
quantities are computed numerically with the LPTR. 

To derive our analogue of the JDFKS fitting formula 
from the ESMR, we use equation (|60|l to convert equation 
(|27|l to b^ . Note that PNG alters the mean ionized fraction 
relative to the Gaussian case at a fixed redshift. We will 
therefore use the super-scripts "NG" and "G" on the mean 
ionized and neutral fractions to distinguish the two cases. 
Taking this caveat into account, equation (|27|) implies 



fNG 



^NG -G 



^H 



£i^G 



h^x- 


'-'min 

6 


51 


-1^ 



'-'minf^a: 



3(5c5n 



51 



y2. 

nun 



+ 25cbj. 



SI 

Almin(fc) 



(63) 



Note that this expression applies to PNG with general bis- 
pectra. In order to compare this result to the JDFKS for- 
mula [equation (|61fl ] we restrict equation (|63|) to the local 
template, where we can make the substitution. 



25,S°^^SLR 



3/nl&3^ 5c 



ilm-n.Q 



Mmin(fc) ' " giO)D{z)k^Tik)- 

Applying this substitution to equation (|63[l yields 



(64) 



fNG 
bx = 


.=NG ;=G r ^ o{3) 


■^min 


^V 


(/ -01 


+ 3/nl6.^(q^^ 


SlmHa y 
J{z)k^T{k) j ■ 



+ 



(65) 

Our analytical result differs from the JDFKS fitting 
formula in several important ways: (i) The second term in 
equation (|65|) is the scale-independent correction to bx due 
to PNG, which JDFKS neglected. For the local template, 
this correction contributes up to a few percent to bx for 
fc ~ 0.1 Mpc~^. As an example, in the ESMR model at 
2 = 10.2 with /nl = 50 (Sp = 0.51, xf = 0.48, b°=4.15), 
we find the fractional contribution of the scale-independent 
term to be 2.5, 1.6 and 0.06 percent at fc = lO"'^, 10~^, 
and 10""^ Mpc~^ respectively, (ii) The critical density for 
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halo collapse, 5r_ ~ 1.686, appears in the last term of (|63p 
(the scale-dependent term) and not 5^, ~ 1 as assumed by 
JDFKS. (iii) Equation (|63|l has an overall factor involving 
mean ionized and neutral fractions which is absent in the 
JDFKS formula, (iv) It is 6^ that sets the amplitude of 
the scale-dependent term in our formula, and not (6^ — 1) 
as assumed by JDFKS. They assumed the latter because 
the non-Gaussian scale-dependent Eulerian halo bias is pro- 
portional to (6 — 1), where h is the Gaussian (Eulerian) 
halo bias. This (6 — 1) form in the halo bias comes from 
the spherical collapse model, in which the Lagrangian and 
Eulerian halo bias parameters are related by & = &l + 1 
(|Mo fc Whit jl 19961 ). However, as we discussed in ^2.3.11 the 
Lagrangian and Eulerian ionized fraction bias parameters 
are equivalent, and since bx = (xi/xH)bx, so are the La- 
grangian and Eulerian h^. The factor of (6^ — 1) should 
therefore not appear in the formula. 

We now explore the numerical differences between equa- 
tions (|61|) and H63p . Figure[Tn]shows the comparison between 
three calculations: (i) The purely numerical LPTRl calcu- 
lation (solid), using equation (|6Up to convert h^§ii to h^'^ , 
(ii) Our ESMR prediction (dashed), equation (|63p . and (iii) 
The JDFKS fitting formula (dot-dashed), equation (|6H) . In 
calculations ii and iii, we insert the numerical values of &^ 
from the LPTRl calculation into equations (|63|) and (|6ip 
respectively (like the "hybrid" calculations described in the 
last section), since our goal here is to determine how well 
these equations capture the mapping between h^ and V^ ■ 
The top panel shows bx as a function of fc at a fixed red- 
shift of z = 10.2, corresponding to a mean ionized fraction 
of « 50%. The bottom panel shows the redshift evolution 
for a fixed k = 0.01 Mpc"\ While the ESMR prediction 
matches the LPTR results well, the JDFKS fitting func- 
tion always yields a significantly lower amplitude of V^ ■ 
The differences can be quantified in terms of /nl, which 
sets the amplitude of the scale-dependent term. As a sim- 
ple, "back-of-the-envelope" illustration, if we assume that 
/nl = 50 is the "true" model, and that the LPTRl curve in 
the top panel of Figure [TO] represents the "measured" V^ at 
Xi ~ 0.5, then an observer using the JDFKS fitting formula 
would falsely infer an /nl ~ 100, assuming that 6^ is known 
exactly. The bottom panel of Figure [TO] shows that the JD- 
FKS formula also does not reproduce the redshift evolution 
of V^ ■ At high redshifts, the JDFKS formula becomes neg- 
ative due to 6^ dropping below unity at those epochs. This 
behavior of the JDFKS fitting formula is not observed in 
our results. 



6 SUMMARY AND DISCUSSION 

We have presented a first-principals investigation on the ef- 
fects of primordial non-Gaussianity on the large-scale struc- 
ture of reionization. We employed two methods that are in- 
dependent in how they model reionization: 1) An extension 
of the analytical excursion-set model of reionization (ESMR) 
oflFurlanctto ct al, (200J) to include PNG, which allowed us 
to derive analytical expressions for the ionized density bias. 
2) The linear pertu rbation theory of reionization (LPTR) of 
IZhang et al.l ((20071), which has the advantage that it directly 
solves the ionization rate and radiative transfer equations, 
allowing us to explore a range of recombination clumping 



factors, and two distinct models of photon-production rates 
in galactic sources. Our main results can be summarized as 
follows: 



• Equations (jlSp and (|19|) inserted into equation (|12|l 
constitute our extension of the analytical ESMR of 
iFurlanetto et al.l ([200J) to include non-Gaussian initial con- 
ditions with general bispectra. 

• In our non-Gaussian extension of the ESMR, equation 
H26I) gives the Gaussian ionized density bias, while equa- 
tions (|32l) and (|33|) give the non-Gaussian scale-independent 
and -dependent ionized density bias corrections. These ex- 
pressions follow from the assumption that, on large scales, 
fluctuations in the ionized density fleld follow fluctuations in 
the source distribution. The ESMR predicts that the ionized 
density bias is scale-independent in models with Gaussian 
initial conditions, while for non-Gaussian initial conditions 
the bias acquires a scale-dependent correction which scales 
as 1/fc^ in the local template for small k. 

• Numerical calculations using the LPTR conflrm that 
the ionized density bias in the small-fc limit (fc < 0.1 Mpc^^) 
is scale-independent in models with Gaussian initial condi- 
tions, and strongly scale-dependent for local PNG, reflect- 
ing the impact of PNG on the source bias. Moreover, the 
imprint of the source bias on the ionized density bias per- 
sists throughout the reionization epoch. We attribute this 
characteristic of our models to the mean free path of UV 
photons through the IGM, which remains small relative to 
the scales on which PNG affects the halo bias, until just 
before reionization ends. 

• While the simple one-parameter ESMR model devel- 
oped in this work cannot capture the detailed evolution 
of even the Gaussian ionized density bias computed in the 
LPTR, we found that equations (|32|l and (|33p provide a rea- 
sonably accurate map between Gaussian and non-Gaussian 
ionized density bias parameters. As shown in Figure [9] these 
formulae work best at earlier times, corresponding to lower 
mean ionized fractions, and in models with lower recombi- 
nation rates. 

• Equation (|63p gives our prediction for the "bias of ion- 
ized" regions deflned by JDFKS, which is related to the ion- 
ized density bias used in this work by equation (|60|l . This 
result can be compared directly to the JDFKS fltting for- 
mula given in equation H61|) . which they used to forecast 
constraints on /nl by future measurements of the EoR 21cm 
power spectrum. We note significant differences in the forms 
of these expressions. In a companion paper, we use our re- 
sult to revisit the topi c of constraining PNG with the EoR 
21cm power spectrum (jMao et al.ll2013l ). 

There are several scenarios one could imagine which 
might require extensions to both our analytical ESMR 
predictions and our LPTR calculations. For example, all of 
our models assumed that atomic cooling halos were the only 
photon-sources during reionization. These were not the first 
halos to form stars, however, and there was likely a modest 
contribution to reionization from minihalos at earlier 
times before feedback mechan isms effectively shut them 
down (see e.g. lAhn et al.ll2012l ). We also did not take into 
account the possible "self-regulation" effects of reionization 
on the star formation rates of lower mass atomic cooling 
halos, below the Jeans flltering scale, M ~ 1O^M0. These 
feedback effects can lead to a more complicated reionization 
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history and possibly a more complicated evolution and scale 
dependence of the ionized density bias (for studies on the 
role of feedback during r eioniz ati on, see e.^. [ Shapiro et alJ 



1994 IChiu fc Ostrikeij I2OO0I: iHaiman fc Holder 
Wvithe fc Loebl bOOsT lOnken fc Miralda-Escudel 



Furlanetto fc Loebl l2005l: iKramer et alJ l2006l: llhev et all 



2013I ). 



20071: I Wvithe fc Morales! [20071: I Wvithe fc Cen[ 



Mesinger fc Diikstral 20081: iKulkarni fc Choudhurv 



2007 



2011 



Alvarez fc Abelll20iafAhn et al.ll2012l : ISobacchi fc Mesinger! 



We have also assumed in the case with no feedback that 
the bias of ionizing sources is equal to the halo bias. This 
assumption is often adopted in the literature, but its accu- 
racy is uncertain. One situation in which this assumption 
might not hold is if mergers play a significa nt role in driving 
star fo rmation during reionization (see e.g. iLamastra et al.l 
((2013) for recent observations at lower redshifts). In this 
case, the source bias might be quite different from the halo 
bias, since the bias of merging halos is different from the gen- 
eral halo population. This difference would likely be reflected 
in the ionized density bias, and it is not clear whether our 
cur rent ESMR expression s would apply in such a scenario 
[see lCohn fc Chanel ([20071 ) for an extension of the Gaussian 
ESMR to include the effects of halo mergers]. Scenarios in 
which the source bias is not simply related to the halo bias 
should be explored to gauge their impact on the ionized 
density bias. 

We have considered a range of reionization models in 
this work, but there is an important characteristic in com- 
mon with all of these models. The mean free path of typical 
photons through the IGM is short relative to the large scales 
of interest, until the end of reionization. This characteristic 
is a consequence of the soft source spectrum we adopted. As 
noted above, this feature is crucial in preserving the imprint 
of the source bias on the large-scale ionized density bias. It is 
therefore important to consider scenarios in which this con- 
dition might not hold. One such scenario is that of a hard 
source spectrum, i.e. a signiflcant contribution from X-ray 
photons, which have a significantly larger mean free path 
relative to UV photons, and may therefore act to suppress 
the amplitude of the ionized density bias on larger scales, 
approaching the scales of interest for the effects of PNG. 
Further numerical work should be devoted towards explor- 
ing more detailed reionization models and their implications 
for the analytical results presented here. 

Finally, we have discussed one observable effect of the 
scale-dependent ionized density bias: a scale-dependent sig- 
nature in the EoR 21cm power spectrum. However, the scale- 
dependent ionized density bias can also have an impact on 
secondary temperature and polarization anisotropics in the 
CMB sourced by the EoR. Extensions of this work may also 
involve observational signatures beyond the redshifted 21cm 
background. 

Note: during the preparation of this manuscript, a pa- 
per appeared on the ArXiv preprint archive on the topic 
of co nstraining PNG with the EoR 2 1cm p ower spectrum 
l[Chongc hitnan 201^'). IChongchitnanI ([20131 ) calculated the 
ionized fraction bias (and not the ionized density bias 
co nsidered in ou r work) based on the reionization model 
of Alvarez et al.l (I2OO6I ). We note that the approach of 
IChongchitnanI ( 20131 ) does not yield an analytical expres- 
sion for the non-Gaussian ionized fraction bias, in con- 



trast to the ESMR approach taken h ere. Shortly after this 
manuscript was submitted, a paper bv lLidz et al.l ([2013h ap- 
peared on the archive which expanded upon the work of 
JDFKS with more-detailed semi-numerical simulations of 
reionization with PNG. They also considered the effects of 
foreground subtraction in their forecasts of constraints on 
/nl from the EoR 21cm power spectrum. 
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The non- Gaussian part comes from 
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coll,0 



dS,r 
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where the "0" sub-script on x again denotes the Sn = and 
Sr = limit. We Taylor expand Sx to first order in Sn 
about Sr = 0. The constant term is zero, while the linear 
term is 
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Taking the limit as 5*^ — >■ yields 



7 (3) 
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where 5^4 = (^min) /Smin deuotcs the skewness of density 
fluctuations smoothed on the Mmin scale. As described in 



i]2.3.1l this is the Lagrangian ionized fraction bias, which 
is equivalent to the Eulerian ionized fraction bias. We may 
therefore apply equation Q to arrive at equation p2|) for 
the scale-independent correction to the ionized density bias. 



APPENDIX A: DERIVATION OF THE 
NON-GAUSSIAN IONIZED DENSITY BIAS IN 
THE ESMR 

In this appendix we present our derivation of the non- 
Gaussian ionized density bias from the ESMR. The first step 
is to compute the ionized fraction bias. We will then use 
equation ^ to convert the ionized fraction bias to the ion- 
ized density bias. We begin with the non-Gaussian ionized 
fraction contrast, equation pijl . which contains both scale- 
independent and -dependent contributions. In what follows, 
we define the following functions for notational convenience: 



1 



1 



3 [ 5111111 — Sr Sc — Sr 

5 c — ScSr 



Sc — {Sc — Sr) coth 



Sr 



(Al) 
(A2) 



Al Scale-independent terms 



The scale-independent contributions can be further divided 
into the Gaussian term and the non-Gaussian correction. 
The first term on the right-hand side of equation (|31l) yields 
the former. Taylor expanding it to first order about Sr = 
and taking the limit as Sr — >■ gives 



^G_ 2 ain/gii^o 

Sc 91n5min ' 
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(A3) 



A2 The scale-dependent term 

Scale-dependent contributions can come from sub-terms in- 
volving (S^^^^Sr) in the last term of equation pif) . There are 
two such terms. 
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Equation 1)31^ is written in coordinate space whereas the 
scale-dependence of the bias is manifested in Fourier space. 
We em ploy a convenient method used bv lDesiacgues et al.l 
(|2011bl ) for the conversion. Adapting their strategy to the 
current task, we take the cross-correlation, (S^Sr), between 
the ionized fraction contrast and the large-scale smoothed 
density contrast. We then rearrange the equations to pick 
off t he Fourier space b ias par ameter. 

iDesiacQues et al.l (|2011bl ) considered the statistics of 
thresholded regions - regions in the initial density field with 
peak height above some given value. They showed that the 
power spectrum of thresholded regions, which in the Press- 
Schechter formalism can be interpreted as the collapsed frac- 
tion power spectrum, can to first order be expressed with a 
simple bias relation. Adopted to our notation, the linear 
bias relation they found is S{coii{k) — fofcoii(fc)'5niiii(fe). Not- 
ing that Sicoii = Sx, and fefcoii = bx in the ESMR, we can 
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write 5x{k) = hx(k)5iain{k). We use this relation to write 
the left-hand side of (1311) as 



{5x5r) 



(27r)3 



ho,{k)Mmiu{k)MR{k)P4,{k). (A9) 



If we can now rewrite the appropriate terms in the cross- 
correlation of 5r with the right-hand side of (|3H) in a similar 
way, we can simply read off the bias coefficient. 

Let us first consider equation (|A7|I . We Taylor expand 
it to first order about 5r, = 0, neglecting the constant term, 
since it will vanish upon taking the cross-correlation with 
5r. Taking the limit Sr — >■ 0, and using the property. 



we can write the cross-correlation of d\ with 5r as 



(AlO) 
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Conveniently, 
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APPENDIX B: DETAILS ON THE SOURCE 
EMISSIVITY IN THE LPTR 

Here we derive expressions for the spatially averaged 
source emissivity and its first-order perturbation, which 
appear in equations (|43|) and (|49|l respectively. For con- 
venience, we define the quantity nu,cou{'r',Muiin,R,'il) = 
nH{r,R,rj)fcou{Mniin,R,SR,rj), representing the density of 
hydrogen collapsed into ACHs, smoothed over scale R. Note 
that TiH.coii appears in equation (|4Up . and its time derivative 
appears in equation (|39]|. 

We write nn.coii in terms of a spatial average and a 
first-order perturbation as follows: 



WH,co1i(t', Af,nm,-R, r?) = Wh 



fcoll,o{Mn-i[n, ri) 



+ All,collir, Mmln, R,v) 



(Bl) 



Referring back to equations (|38|l . (|39|) . and (|40|l . the 
quantity ^s which appears in equation (|43|l can sim- 
ply be read off from equation (|B1|) : in source-model A, 
^f = 'y^ixv)d[fcoii,o{Mn^in,ri)] /dri, while in source-model 
B, ^s = 7''(a:-)/coii,o(M„in,77). 

If we write the smoothed hydrogen number density as 
nnir, RjTj) — fT,H(77) [1 -I- 5_B(r, ?;)], and the conditional col- 
lapsed fraction as /coii(Mmm,-R, ^ii,??) = /coii,o(Mmin,77) -I- 
A/coii(Af,nin,-R, (5fl,?7), then to leading order Ah,co11 = 
5R{r,r])fcou,o{Mn^in,ri) + A{coiiiMinin,R, Sr,!]), aud its 
Fourier transform is 



(diSR,) 






(A13) 



In a similar manner it is straightforward to show that 
{d25R) = 0. 

The mixed correlator, (iJ^in^fl), can be written in a form 
similar to equation (IA9II [see disc ussion leading up to equa- 
tion (9) of lD'Aloisio et all (|2012l )]. 



(^min'^fl) = 






-MR{k,z)P4,{k)4:Smln J^'^Uk). (A14) 



(3) 



Plugging this into equation (|A13[) . in combination with 
equation (|A9|) . yields 



A6i'')(fc) = 25,6^-^-'"^''^ 



Xmin(fc) ' 



(A15) 



where the form factor .7^min(^) i^ defined in equation (|34p. 
W e note that this expressio n is equivalent to equation (43) 
of iDesiacgues et al.l (|2011br ) for A'^ = 3. Their linear bias 
of thresholded regions is equivalent to our collapsed fraction 
bias, defined in equation (|lip . Since the ionized fraction bias 
is equivalent to the collapsed fraction bias in the ESMR, the 
former is also equivalent to the lin ear bias of thresholded 
regions in IDesiacgues et al.l (|2011a ). hence the correspon- 
dence between their equation (43) and equation (|A15|l . Fi- 
nally, equations Q and (|A15|) yield equation (|33|l for the 
scale-dependent correction to the ionized density bias. 



AH,coii(k, ry) = /coii,o(77)<5(k,77) + Afcoii(k,r7), 



(B2) 



where from here on we suppress the Mmin dependence for 
brevity (it should be understood that the collapsed fraction 
refers to the fraction of mass in halos above the ACH thresh- 
old), and the dependence on R has been dropped under the 
assumption of small k. The quantity Afcoii may be written 
in terms of the collapsed fraction bias, 

Afcoii(k,77) = fefcoii(k,77)(5(k,?7)/coii,o('7), (B3) 

from which equation HB2|) becomes 



i(k,J?) = fcoiifiiv) 



l + bfcoii(k, r;) 



S{KV)- (B4) 



The quantity As appearing in equation (|49p can now be read 
off, yielding equations (|53|) and (|54|) in source-models A and 
B respectively. 

The collapsed fraction bias parameters which enter 
these expressions follow trivially from the results of ap- 
pendix |X] since, in the ESMR, the ionized fraction bias is 
assumed equal to the collapsed fraction bias. The Gaus- 
sian collapsed fraction bias is therefore equation (|A3|I . while 
the scale-independent and -dependent terms from PNG are 
()A6|I and (|A15|I respectively, each with the replacement 
bx -^ bfcoii- 
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